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ABSTRACT 

Many young stars reside within the central half-parsec from SgrA*, the supermassive black hole in the Galactic 
Center. The origin of these stars remains a puzzle. Recently, Hansen and Milosavljevic (2003, HM) have argued 
that an Intermediate-Mass Black Hole (IMBH) could have delivered the young stars to the immediate vicinity of 
SgrA*. Here we focus on the final stages of the HM scenario. Namely, we integrate numerically the orbits of 
stars which are initially bound to the IMBH, but are stripped from it by the tidal field of SgrA*. Our numerical 
algorithm is a symplectic integrator designed specifically for the problem at hand; however, we have checked our 
results with SYMBA, a version of the widely available SWIFT code. We find that the distribution of the post- 
inspiral orbital parameters is sensitive to the eccentricity of the inspiraling IMBH. If the IMBH is on a circular 
orbit, then the inclinations of numerically computed orbits relative to the inspiral plane are almost always smaller 
than 10 degrees, and therefore (a) the simulations are in good agreement with the observed motions of stars in 
a clockwise-moving stellar disc, (b) the simulations never reproduce the orbits of stars outside this disc, which 
include those in the second thick ring of stars and the randomly oriented unrelaxed orbits of some of the S-stars. 
If the IMBH's orbital eccentricity is e = 0.6, then approximately half of the stars end up with orbital inclinations 
below 10 degrees, and another half have inclinations any where between and 180 degrees; this is somewhat closer 
to what's observed. We also show that if IRS 13 cluster is bound by an IMBH, as has been argued by Maillard 
et. al. 2004, then the same IMBH could not have delivered all of the young stars to their present location. 

Subject headings: black holes, massive stars, orbits 



1. introduction 

A few tens of bright, young, and massive stars have been ob- 
served in close proximity (0.001 — 0.5 pc) to the radio source 
SgrA*, which is associated with the supermassive (4 x 1O 6 M ) 
black hole in the Galactic Center 1 (see Genzel et. al. 2000, Ghez 
et. al. ,2003, Shodel et. al. 2002,2003 for recent work). The 
apparent youth of these stars implies two logical possibilities 
for their birth history: (1) either they were born in situ, in the 
immediate vicinity of SgrA*, or (2) they were born some dis- 
tance away and then delivered to SgrA* within a few million 
years after their birth. In the first scenario, the in situ birth 
would occur in a strong tidal field of SgrA*, and the density of 
the star-forming gas would have to be several orders of mag- 
nitude greater than that anywhere else in the Galaxy (Phin- 
ney 1989, Sanders 1992, Morris 1993). The most likely ge- 
ometry for the star-forming gas is a dense disc, which is ei- 
ther compressed by the central explosion powered by SgrA* 
(Morris, Ghez, and Becklin 1999), or, more plausibly, which is 
massive enough to become self-gravitating and fragment into 
stars 2 (Levin and Beloborodov 2003, Nayakshin, Cuadra, and 
Sunyaev 2004, Milosavljevic and Loeb 2004, Nayakshin and 
Cuadra 2004). The in situ formation is not the focus of this 
paper, although we will briefly return to it in our concluding 
remarks. 

The second scenario, which is the focus of this paper, calls 
for rapid delivery of the young stars to the supermassive Black 
Hole. The initial version of it was proposed by Gerhard 2001, 
and investigated in more detail by McMillan and Portegies 
Zwart 2003, and Kim and Morris 2002. Gerhard's basic idea 
was that if a massive (~ 1O 6 M ) cluster of stars was born 



10 — 30pc away from SgrA*, then the dynamical friction in 
the central bulge would bring the cluster towards SgrA* on the 
timescale of a few million years. During the inward migration, 
the massive stars in the cluster would sink to the cluster's center 
and form a compact core. Gerhard has argued that even though 
the envelope of the cluster gets stripped by the tidal field a few 
parsecs away from SgrA*, the dense cluster core would only 
be tidally disrupted within the central parsec. However, even 
that is not close enough: the observed stars reside 0.001 — 0.5 
pc from SgrA* and the core density would have to increase by 
orders of magnitude before it could get so close to SgrA*. 

Recently, Hansen and Milosavljevic 2003 (from now on, 
HM) have suggested a fix to Gerhard's scenario which would 
enable the stars to be delivered to the observed proximity from 
SgrA*. They have argued that if an Intermediate-Mass Black 
Hole (IMBH), with the mass of 10 3 — 10 4 M Q , was positioned 
at the center of the cluster, than it would provide some extra 
gravitational binding for the cluster core. This scenario is sup- 
ported by an observation of the seemingly bound mini-cluster 
IRS 13 at ~ O.lpc from SgrA* (Maillard et. al. 2004, from here 
on MPSR). MPSR have argued that the presence of an IMBH 
at the center of IRS 13 is the most plausible explanation of 
how the mini-cluster keeps itself together in the strong tidal 
field of SgrA*. On the theoretical side, a number of investi- 
gations argue that IMBH would plausibly form in a dense core 
of a young very massive cluster (Portegies Zwart et. al. 2004, 
Gurkan et. al. 2004). The cluster+IMBH inspiral would consist 
of two stages: the first one, in which most of the most of the 
cluster is tidally disrupted at about a parsec from SgrA* and 
the IMBH is released together with a few tens of massive stars 
tightly bound to it, and the second one, in which the IMBH con- 



'For convenience, we shall refer to the supermassive black hole simply as SgrA*. 

2 We note that it is by no means certain that the disc would circularize before forming stars (J. Goodman, private communications). 
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tinues its inspiral while the stars are peeled off one-by-one by 
the SgrA* tidal pull, and are eventually put into nearly Keplerial 
orbits around SgrA*. The first stage was recently investigated 
numerically by Kim and Morris 2003 (KM) and by Gurkan and 
Rasio 2004. KM have found that very few stars survive bound 
to IMBH once it gets to the central parsec. However, one of the 
numerical limitations in KM (pointed out by the authors) was 
the large softening length for gravitational interactions between 
the stars. Thus KM did not model faithfully the internal dynam- 
ics of the cluster. A more exact Monte-Carlo-based simulation 
of Gurkan and Rasio 2004 has found that the number of stars 
bound to IMBH at the end of the first stage is consistent with 
the number of young stars in the central half -parsec. 

Here we focus on the second stage — the final inspiral of the 
IMBH. We investigate numerically the stellar orbits produced 
by such an inspiral, and compare them to observations. In the 
next section we describe the algorithm which was used for the 
numerical orbital integrations; this section may be skipped by a 
reader uninterested in technical details 3 . In section 3 we present 
our numerical results, and in section 4 compare them to the 
data. We find that while the inspiral on a circular orbit explains 
nicely the clockwise-moving stellar disc in the Galactic Cen- 
ter (Levin and Beloborodov 2003, Genzel et. al. 2003), it does 
not reproduce the thick counterclockwise stellar disc and the 
randomly oriented tight orbits of the S-stars (Ghez et. al. 2003, 
Shodel et. al. 2003). The inspiral on an eccentric orbit does 
push some orbits out of the inspiral plane, but is still not en- 
tirely consistent with the current data. We also show that the 
proposed IMBH at the center of IRS 13 almost certainly could 
not have delivered all of the young stars to the Galactic Center. 
Finally, we speculate how the HM scenario may be modified to 
account more closely for all of the kinematic data. 

2. OUR NUMERICAL ALGORITHM: SYMPLECTIC INTEGRATOR IN 
THE EXTENDED PHASE SPACE 

The stars in our system will experience repeated close en- 
counters with the IMBH even after they get tidally stripped 
from it, and we expect the resulting orbits to be very chaotic 
and sensitive to the initial conditions. Therefore, our study will 
be statistical in nature: we are seeking to find the distribution 
of orbits which are produced by the IMBH inspiral. Symplectic 
algorithms (SA; see Yoshida 1990) are particularly suitable for 
statistical studies, since they generate trajectories which con- 
serve the phase-space volume exactly. Integrators which use 
SAs have excellent long-term behavior, since the Hamiltonian 
mapping used in a SA is integrated exactly and is only slightly 
different from a true Hamiltonian mapping of the system. How- 
ever, for our purposes the simplest symplectic codes have a dis- 
advantage since they do not handle well the close encounters 
between a star and IMBH (see Preto and Tremaine 1999 for dis- 
cussion). The origin of this limitation is the reduced timestep 
which is used to integrate through the encounter; reduction of 
the timestep leads to change of the SA Hamiltonian and thus 
may cause spurious energy changes. To overcome this prob- 
lem, Mikkola (1997), Preto and Tremaine (1999), Mikkola and 
Tanikawa (1999), and Mikkola and Wiegart (2002) have devel- 
oped extended phase-space Symplectic Integrators which use 
fixed timestep but in an extended phase space. The basic idea 
for extending the phase space of the Hamiltonian system goes 
back to Poincare, and can be briefly described as follows: 

Let H(q,p,t) be a time-dependent Hamiltonian, where q and 

3 In our view, however, this algorithm is of interest in its own right and may be 



p stand for spatial coordinates and their conjugate momenta re- 
spectively and t is the time. Let us extend the phase space of 
the system by another coordinate qo and its conjugate momenta 
Pq. Consider now a new Hamiltonian T in the extended phase 
space (Mikkola 1997): 

r(q,q ,p,p ) = g(q,q ,p,p ) [H(q,p,q ) + p ] , (1) 

where g is some smooth differentiable function. This Hamil- 
tonian depends only on the coordinates and momenta of the 
extended phasespace; thus any trajectory in the extended phase 
space which satisfies Hamiltonian equations will conserve T. 
Therefore, if a particle begins its motion on the hypersur- 
face po = -H(q,p,qu), it will always stay on this hypersurface. 
Moreover, the trajectory of this particle [q(T),p(T),qo(T)] will 
trace the values [q,p,t] corresponding to the trajectory of the 
original Hamiltonian system; here r is the time variable which 
marks evolution in the extended phase space. The two time 
variables are related by 

dt = g(q,q Q ,p,-H)dr. (2) 

The idea then is to integrate the system numerically in the ex- 
tended phase space using the constant timestep dr, but choos- 
ing g so that the real timestep dt is small during a close en- 
counter. 

To integrate the motion in a symplectic way, we need the 
Hamiltonian T to be separable into two parts, 

r = r A +r B , (3) 

so that the motion due to each of Fa and Lb would be integrable 
analytically. Once such separation is found, the system can then 
be integrated via a generalized leapfrog step: 

U(dr) = U A {dT/2)U B (dT)U A (dT/2), (4) 

where UA(dr) and Ub(At) represent the forward motion in the 
extended phasespace by the time interval dr under the action of 
L4 and Tb, respectively. Preto and Tremaine 1999 and Mikkola 
and Tanikawa 1999 have independently found the form of the 
generalized Hamiltonian which is the sum of two analytically 
integrable parts: 

T = f(p +K)-f(-U[q, qo ]), (5) 

where K and U are the kinetic and potential energies, respec- 
tively, and / is an arbitrary smooth function. The Hamil- 
tonian in Eq. is obviously of the form in Eq. Q, with 
g = T/(po+K + U). When po = -H, i.e. at the hypersurface 
representing the original Hamiltonian system, g = /'(—{/). The 
choice fix) = logx possesses the following surprising and useful 
property: the generalized leapfrog step for a two-body problem 
follows exactly the correct conic-section orbit, albeit with some 
time error. The real time and the time in the extended phase 
space are then related by 

dt = dr/(-U). (6) 

We can now see the advantage of using the extended phase 
space. In the extended phase space, we can integrate the motion 
by using the generalized leapfrog in Eq. (0} with the timestep 
fixed, which is required for energy conservation. However, the 
timestep in real phase space will not be constant: when the 

i in other contexts. 
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potential energy becomes large during a close encounter, the 
timestep in real space becomes small and the precision of real- 
space integration increases. 

For our problem however, the generalized Hamiltonian in 
Eq. (jSJl is not optimal. We are interested in a restricted 3- 
body problem where a star is interacting with two much heavier 
Black Holes. The mass ratio between SgrA* and the IMBH is 
1000 or more 4 and our star is initially bound to IMBH. The 
potential energy of the star per its unit mass is given by 



U = - 



M 



■n 



(7) 



where r is the position vector of the star, M and f\ are the mass 
and the position vector of SgrA*, and m and 7i are those of the 
IMBH. In our simulations, the star is treated as a massless parti- 
cle and the orbits of the two black holes are introduced by hand, 
with the inspiral prescription motivated by analytical calcula- 
tions. We see from Eq. Q that the potential energy of the star 
is dominated by its interaction with SgrA*, even when it is well 
inside the Roche Lobe of the IMBH. Thus if the Hamiltonian 
of Eq. (jSJi is used to generate the integrator, the following para- 
doxical situation may occur: the star may be on a bound orbit 
around the IMBH or it may be experiencing a close encounter 
with the IMBH, yet the choice of timestep will be dominated 
by the star's interaction with SgrA* and not with the IMBH. 
Clearly, to make the integrator efficient, we need to reduce the 
relative contribution to the timestep from the star-SgrA* in- 
teraction. The following generalized Hamiltonian, found by 
Mikkola and Wiegert in 2002, does the trick: 
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Here all the position vectors are measured relative to the 
barycenter of black-hole binary, which is very close to SgrA*, 
a -c 1 is a positive number, and the term 5 is given by 



5 = Mx 



1 



l 



(9) 



Since we treat the star as a test particle and introduce the or- 
bits of both black holes by hand, the vectors 7\ and 7i are 
time-dependent and therefore are prescribed functions of qo. 
Both components of the Hamiltonian in Eq. are easily in- 
tegrable: the first term on the RHS describes Keplerian motion 
around the barycenter and the second term describes the motion 
where all positions stay constant and all momenta change lin- 
early with time; therefore we can efficiently construct the gen- 
eralized leapfrog operator for this Hamiltonian. The timestep 
in real phase space is then given by 



dt = dr x 



-I 2 



(10) 



If one chooses a = (m/M) 2 / 3 , then at the IMBH Roche sur- 
face the star-SgrA* interaction and star-IMBH interaction con- 
tribute approximately equally to the timestep, and the contribu- 
tion from the <5-term is relatively small. 



Mikkola and Wiegert (2002) have proposed to use the gener- 
alized Hamiltonian in Eq. for simulating the motion of near- 
Earth and near- Jupiter asteroids. They have cautioned however, 
that if the asteroid gets too close to the sun, the <5-term in Eq. (|8} 
may become large, and the Hamiltonian may become singular. 
Our simulations have confirmed this. We have found that while 
for circular inspiral the Hamiltonian in Eq. (|8j always works 
well, for an eccentric inspiral some of the simulation runs crash 
because the Hamiltonian becomes singular. However, we have 
identified a way to fix this singularity problem by a slight mod- 
ification of both parts of the generalized Hamiltonian: 

1- 



(11) 



r = io g PQ+K-—— + —T 



2r 2 



log 



■'"2 1 



The above Hamiltonian is of the form in Eq. Q and, moreover, 
the two terms on the RHS are still separately integrable. The 
latter statement is not trivial for the first term on the RHS; how- 
ever, recall that for a particle motion in a spherically symmet- 
ric potential V(r) the effective potential for the radial motion 
is given by V(r) + (l/2)L 2 /r 2 , where L is the particle's angular 
momentum. Thus adding the term (l/2)e 2 /r 2 to the Hamilto- 
nian is equivalent to rescaling of the particle's angular momen- 
tum in the radial equation of motion: L — > y / L 2 + e 2 , and solving 
the equations of motion with this term can be reduced to solv- 
ing a purely Keplerian problem 5 . It is straightforward to show 
that by setting the value of e so that e 2 /2 > mr™ one can avoid 
the singularity in the generalized Hamiltonian; here r^ax is the 
maximal distance from the IMBH to the barycenter during the 
simulation run. 

In this work, we have used the Hamiltonians in Eq. (|8j and 
Eq. (II 1 fc to simulate the circular and eccentric IMBH inspirals, 
respectively. We have also, as a check, have run test simulations 
with SYMBA, a version of widely available SWIFT code with 
the added ability to resolve close encounters between massive 
objects (Duncan, Levison, and Lee 1998). We have confirmed 
that our integrator and SYMBA give results which are consis- 
tent with each other. The efficiency of our code has allowed 
us to run thousands of inspiral simulations and obtain a distri- 
bution of the orbital parameters of the stars at the end of the 
inspiral. In the next section, we discuss the results of our simu- 
lations. 



3. RESULTS 



3.1. 



Parameters of the inspiral. 

We show in the Appendix that the eccentricity of the inspi- 
raling IMBH remains nearly constant for the observed density 
profile of the stellar cluster. The semimajor axis evolves with 
time as 

a(f) = floexp(-f/f ), 
where to is the characteristic inspiral time given by 



(12) 



M 1 



87rVGm / 9r L5 log(A) 



(13) 



where M is the mass of SgrA*, m is the mass of the IMBH, p is 
the mass density of stars in the central cluster, and ln(A) is the 



4 The mass of the IMBH is limited from above by the mass of massive stars participating in the core collapse of the IMBH's parent cluster, about 0.1% of the 
cluster's mass. 

5 We will discuss elsewhere how we have implemented this procedure in practice. 
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Coulomb logarithm. For circular orbit, this can be rewritten as 



'orb 



167T 2 l0g(A) 
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(14) 



where r orb is the orbital period of the IMBH. Genzel 
et. al. (2003) have measured the density profile of the central 
cusp to be 

p~ lO 6 (r/loV' 3 M pc- 3 . (15) 



For a realistic IMBH M/m ~ 1000. Therefore for the stars 
located a few arcseconds from SgrA* the relevant inspiral 
timescale is between a few hundred and a few thousand or- 
bital periods. In our simulations we use fiducial values of the 
mass ratio M/m = 1000 and the inspiral timescale of one thou- 
sand initial orbital periods of the IMBH. Our inspiral timescale 
remains constant over the simulation run; this approximation 
would be exact if the cusp density profile scaled as p oc r~ l 5 (see 
e.g. Appendix of Gould and Quillen, 2003). We spot-check our 
results by varying the inspiral timescale and making sure that 
our qualitative conclusions do not change. 

3.2. Circular inspiral. 

Figure 1 shows a pictorial representation of a typical circular 
inspiral. The points (x,y) on the figure represent ten thousand 
snapshots of the x and y coordinate of the star. The orthogonal 
rotating axis x and y belong to the inspiral plane and are chosen 
so that the SgraA* and IMBH are always on the x axis. From 
Figure 1 we can trace the dynamical history of the star as the 
inspiral proceeds. Initially the star is bound to the IMBH which 
is climbing towards the center on the x-axis. Then the tidal field 
of SgrA* disrupts the star-IMBH binary, the star escapes from 
the IMBH via the inner Lagrange point and spends some time 
inside of the IMBH orbit. Then as a result of a close encounter, 
the star is flung outside of the IMBH orbit. The star then con- 
tinues to have encounters with the IMBH, and its eccentricity 
and inclination get kicks during the encounters. Eventually, the 
orbit of the IMBH is shrunk enough so that the encounters do 
not occur any more and the star's eccentricity and inclination 
stay virtually constant. Figure 2 and 3 show how the star's ec- 
centricity and inclination evolve with time. 

Using CITA's Mackenzie cluster, we have performed one 
thousand simulation runs, such as the one shown in Figures 
1,2, and 3. The initial conditions were different for each run; 
we have chosen the stars in all the runs to have the same ini- 
tial Jacoby integral. In Figure 4 we show in the (e,i) plane the 
scatter plot of the final eccentricities and inclinations. Because 
of multiple close encounters, the orbits are highly chaotic: even 
a slight change in the initial conditions or in the timestep of 
the integrator leads, after some time, to a very different tra- 
jectory and thus to a different final eccentricity and inclination 
relative to the inspiral plane. One subtle issue is the choice of 
timestep 6 . Usually, one can make sure that the chosen timestep 
is appropriate by running the simulation with the fraction of the 
chosen timestep and checking if the simulation results stays the 
same. In our case however, the system is chaotic and the in- 
dividual trajectory will always be different after some time if 
a different integration timestep is chosen. However, it is not 
the individual trajectories that we are after; rather, we want to 
know what distribution of the orbital parameters one should ex- 
pect in the inspiral scenario. Therefore, we have repeated the 
1000 runs with the same initial conditions but with the timestep 

6 Here we mean the timestep dr in the extended phase space; see the previous 



of 0.3 of the original timestep, and checked if the two sets of 
(e, i) points could belong to the same two-dimensional distribu- 
tion. The first thing to check was whether the averages agreed; 
they did. However, we wanted a more complete test. For one- 
dimensional data, Kolmogorov-Smirnov test is both popular 
and mathematically well justified. For two-dimensional data, 
there exist only semi-empirical algorithms based on Monte- 
Carlo simulations: see Peacock (1983) and Fasano and Frances- 
chini (1987). We used the procedure outlined in the latter to 
check whether the two (e, z) data sets were compatible with the 
same distribution; they were. 

The remarkable thing about the scatter plot in Fig. 4 is that 
while the stars, on average, end up with significant eccentric- 
ity, more than 99% of them have inclinations smaller than 10 
degrees. We make another 1000 runs with the inspiral rate re- 
duced by a factor of 10, we find that while the average eccen- 
tricity has increased significantly, and the distribution of incli- 
nations remains similar, with only a few greater than 10 de- 
grees. 

3.3. Eccentric inspiral. 

The final inclinations are no longer small when the inspi- 
raling IMBH is on an eccentric orbit. Qualitatively, eccentric 
stellar orbits experience secular torque if the semimajor axes of 
the star and IMBH are misaligned. Hence angular momentum, 
inclination, and eccentricity of the star can experience substan- 
tial secular changes during the course of the inspiral. In Figure 
5 we plot a fraction of stars with inclinations greater than 10 
degrees for a few values of IMBH eccentricity. Each point of 
the graph is obtained by simulating 100 stellar orbits, and the 
stars which become unbound from SgrA* are removed from 
the data set. In Figure 6, we show average inclination of stellar 
orbits after the inspiral, for the same values of IMBH eccentric- 
ity. In Figure 7, we show the scatter plot of final inclinations 
and eccentricities of stars for the case when the IMBH eccen- 
tricity is 0.6; the scatter plot contains 1000 points. We observe 
two groups of unejected stars which are roughly equal in size. 
In one group, the inclinations are less than 10 degrees, and in 
another the inclinations are spread evenly between 10 and 180 
degrees. In Figures 8 and 9 we show the time evolution of in- 
clination and eccentricity for a typical star which ends up with 
high inclination. 

4. DISCUSSION AND COMPARISON WITH THE DATA. 

Current kinematic data indicates that there are two nearly or- 
thogonal stellar discs in the Galactic Center. One of the disks 
consists of clockwise-moving stars at ~ 0.1 pc from the Black 
Hole, and has dispersion in inclinations no larger than 10 de- 
grees (Levin and Beloborodov 2003, Genzel et. al. 2003). Some 
of the stars in this clockwise disc have significant eccentricities 
(Beloborodov and Levin, 2005). The second disc consists of 
stars which move counter-clockwise (Genzel et. al. 2003) and 
its thickness has not been estimated yet; however it is probably 
thicker than the clockwise disc since the planar fit to its stellar 
velocity vectors has a large \ 2 of 3.5. The mini-cluster IRS 13 
with a putative IMBH belongs to this second, thicker disc. 

In our simulations the circular inspiral nicely reproduces the 
kinematics of the thin clockwise disc: dispersion of inclina- 
tions is less than 10 degrees yet some of the stars possess 
significant eccentricities. However, the nearly circular inspi- 
ral would never produce the second population of stars which 



section. 



5 



move counter-clockwise. Eccentric inspiral seems to be more 
promising: while significant fraction of stars remain close to 
the inspiral plain, the rest of the stars get pushed out of the in- 
spiral plain and end up with large inclinations; see Figures 5 
and 7. However, the stars with large inclinations do not have 
a preferred orbital orientation, thus it is hard to imagine them 
assembling a second disc. Upcoming observations will better 
constrain the thickness of the second disc and then it will be- 
come possible to make a more quantitative assessment of how 
likely it is to produce the observed counter-clockwise trajec- 
tories by an IMBH inspiraling through the plane of the clock- 
wise disc. It seems certain, though, that an IMBH in a counter- 
clockwise moving IRS 13 would not be able to produce the thin 
clockwise disc. Thus the presence of IMBH in IRS 13 could 
not account for kinematics for all of the young stars at ~ 0.1 pc 
from SgrA*. 

The IMBH in IRS 13 would also have difficulty generating 
highly compact eccentric orbits of the most central stars like 
SO-2 or SO-16 (see Ghez et. al. 2005). These stars' apoapses 
are an order of magnitude closer to SgrA* than IRS 13. Such 
disparity of stellar and IMBH orbits is never observed in our 
simulations. However, we note from Figure 9 that in the case 
of an eccentric inspiral, the test bodies with high inclination 
undergo secular evolution in which their eccentricity reaches 
values close to 1, and their distance of closest approach to the 
central Hole is much smaller than their semimajor axis. Thus, if 
such a test body were a binary, this binary could get disrupted 
by SgrA*, and one of the components of the binary could re- 
main tightly bound to SgrA*. The binary-disruption scenario 
for formation of the SO-stars was already proposed by Gould 
and Quillen (2003) in a different context. 

We thank Scott Tremaine for introducing us to symplectic 
integrators in the extended phase space, and Carlo Contaldi for 
showing us how to use the CITA McKenzie cluster. We have 
benefited from discussions with Andrei Beloborodov, Jeremy 
Goodman, Brad Hansen, Chris Matzner, Milos Milosavljevic, 
and especially with Norm Murray. Our research was supported 
by NSERC. 

APPENDIX: IMBH ECCENTRICITY EVOLUTION DURING THE 
INSPIRAL. 

The orbit of an IMBH in the Galactic Center shrinks due to 
dynamical friction, and its eccentricity may also evolve with 
time. The details of this evolution depend on the structure of 
the stellar cusp surrounding the Black Hole. It is instructive 
and convenient to model the cluster by a collection of stars of 
mass with an isotropic distribution function 



f(r lV ) = (3E\ 



(16) 



where E = GM / r-v 2 /2 is the binding energy of the star, and 
and 7 are real numbers. For a Keplerian orbit, the angular 
momentum and binding energy are given byL= y/GMa(l-e 2 ) 
and E = Q.5GM /a, respectively; here e is the eccentricity of the 



orbit. From these relations it is straightforward to derive the 
evolution equation for IMBH eccentricity: 
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(17) 



Here the symbol () represents averaging over the orbital period. 
The acceleration from dynamical friction acting on the IMBH 
with velocity vector v is given by the Chandrasekhar formula 



where 



a = —kv, 



47rG 2 m*mlogA f v , 
k = / dux4iruj{r,u). 



Therefore, 



and 



(dE/dt) = (kv 2 ), 



(dL/dt)=-(k)L; 



(18) 



(19) 



(20) 



(21) 



cf. Eqs (A5) and (A6) of Gould and Quillen (2003). It is conve- 
nient to perform orbital averaging by expressing all the quanti- 
ties through the mean anomaly <f>: 



r = a(\ -<?cos</>), 
, GM l+ecos</> 
a 1— ecos0 

dt/T 0lh = (l-ecos(/>)^. 

2tt 



(22) 



After some algebra, one gets 

dloge l-e 2 
a! log a 2e 2 

where 



(ki/k 2 -l), 



(23) 



1 ,-271 

.2, 



ki = I x L dx I d - (1 +e cos cj})x 2 r (l-e cos <j}) 1 ^, (24) 
Jo Jo 

and 

f 1 , f 27T , 
k2= x dx d4>[l-(l+ecos(f))x ] 7 (l-ecos</>) 7 (l+ecos</>). 

(25) 

We evaluate the integrals numerically and in Fig. 10 we plot 
-dloge/dloga as a function of eccentricity, for 7 = -0.2 and 
7 = 0.25. The two values of 7 correspond to the observed 
density profile p cx r~ L3 and the Bahcall-Wolf density profile 
p oc r" L75 . We see from the graphs that in both cases eccentric- 
ity changes by a factor less than 2 as the semimajor axis shrinks 
by 3 orders of magnitude. Hence, in our simulations we are 
justified in keeping the eccentricity of IMBH orbit constant. 
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FIG. 1 . — 10000 snapshot of the steallar position in the IMBH inspiral plane. Rotating coordinates axes are chosen so that both IMBH and SgrA* are on the x-axis, 
and their origin coinsides with the baricenter. The star is originally bound to the IMBH; this corrensponds to the dense dark region of the scatter plot. The star 
escapes through the inner Lagrange point, then gets flung outside by a close encounter with the IMBH. A few close encounters follow, and then the orbit stabilises 
once the black-hole binary is shrunk away from the stellar periapse. Units on x and y axes are arbitrary. 



8 



2000 



4000 



o 
o 
o 







6000 
3 



2000 



4000 



6000 



q0 



FIG. 2. — 10000 snapshots of eccentricity evolution for the star in Figure 1. The eccentricity is evaluated relative to the baricenter. When the star is bound to the 
IMBH and comes close to it, the eccentricity relative to the baricenter may become greater than 1. Time units: initial IMBH orbital period is 2n. 
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FIG. 3. — 10000 snapshots of inclination evolution for the star in Figure 1. The inclination is evaluated relative to the baricenter. When the star is bound to the 
IMBH and comes close to it, the inclination relative to the baricenter may become large. Time units: initial IMBH orbital period is 2-7T. 
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Fig. 4.— 



Scatter plot of eccentricities and inclinations for a thousand stars after a circular inspiral. 



11 




0.8 



0.6 



0.4 



0.2 



0.4 

IMBH eccentricity 



FIG. 5. — A fraction of stars with inclination > 10 degrees plotted vs IMBH eccentricity. Each point is obtained from 100 runs. 
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Fig. 6.— 



Average inclination plotted vs IMBH eccentricity. Each point is obtained from 100 runs. 
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FIG. 7. — Scatter plot of eccentricities and inclinations for a thousand stars after an inspiral with IMBH eccentricity of 0.6. 
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FIG. 8. — Inclination evolution of a typical high-inclination star; IMBH eccentricity is 0.6 



15 



2xi0 4 4xl0 4 6xl0 4 

T 1 1 1 1 1 1 1 1 1 1 



J I I _J I I I I I I I 

2xl0 4 4xl0 4 6xl0 4 



qO 



FIG. 9. — Eccentricity evolution of the same high-inclination star as in Figure 8. 




FIG. 10. — Plot of the IMBH eccentricity growth rate dloge/d\og(l/a) vs e. Upper curve is computed for the observed density profile of luminous matter in the 
Galactic Center, p oc r~ L3 , and lower curve is computed for a relaxed Bahcall-Wolf density profile, p oc r -1 Js . 



